Efectos del muestreo y aliasing¶

Efectos del muestreo¶

¿Qué le ocurre al espectro de una señal continua cuando la muestreamos?¶

Sin pérdida de generalidad, consideremos para este ejemplo la señal gaussiana

\[ s(t) = \exp \left ( -\frac{t^2}{2\sigma^2} \right) \]

Muestrear es equivalente a multiplicar una señal continua por un tren de impulsos, también conocido como «peineta de Dirac»

\[ \upuparrows(t) = \sum_{m=-\infty}^\infty \delta[t - m / F_s], \]

donde $F_s$ es la frecuencia de muestreo, es decir el inverso entre la separación de los dientes de la peineta

t = np.arange(-5, 5, step=0.0001)
s = lambda t, sigma=0.5 : np.exp(-0.5*t**2/sigma**2)

Fs = Slider(start=0.2, end=5, value=1, step=.01, title="Frecuencia de muestreo")

t_dirac = np.arange(-5, 5, step=1/Fs.value)
create_figure = lambda title : Figure(plot_width=300, plot_height=280, title=title, 
                                      toolbar_location="below", x_range=(-5, 5))
p1 = create_figure('Señal Gaussiana')
p2 = create_figure('Peineta de Dirac')
p3 = create_figure('Señal muestreada')
p1.line(t, s(t), line_width=3)
for p in [p1, p2, p3]:
    p.xaxis[0].axis_label = 'Tiempo [s]'
source = ColumnDataSource(data=dict(t=t_dirac, 
                                    tops=np.ones_like(t_dirac),
                                    sd=s(t_dirac)
                                   ))

p2.vbar(x='t', top='tops', source=source, width=0.05)
p3.scatter('t', 'sd', source=source)


callback = CustomJS(args=dict(source=source, Fs=Fs), code="""
    var t_dirac = [];
    var tops = [];
    var sd = []
    for (let i = -5.0; i <= 5.0; i+=1/Fs.value) {
        t_dirac.push(i);
        tops.push(1);
        sd.push(Math.exp(-0.5*Math.pow(i, 2)/Math.pow(0.5, 2)));
    }
    source.data['t'] = t_dirac;
    source.data['tops'] = tops;
    source.data['sd'] = sd;    
    source.change.emit();
""")


Fs.js_on_change('value', callback)
show(column(Fs, row(p1, p2, p3)))

La transformada de Fourier del tren de impulsos es

\[ \mathbb{FT}[\upuparrows(t)] = F_s \sum_{m=-\infty}^\infty \delta(f - m F_s) \]

es decir otro tren de impulsos pero en frecuencia

Pueden revisar la demostración de la transformada anterior aquí

La señal muestreada es igual a $s(t) \cdot \upuparrows(t)$

Para obtener el espectro de la señal muestreada recordemos que multiplicar dos señales en el dominio del tiempo corresponde a convolucionar sus transformadas de Fourier en frecuencia

\[\begin{split} \begin{align} \mathbb{FT}[s(t) \cdot \upuparrows(t)] &= \mathbb{FT}[s(t)] * \mathbb{FT}[\upuparrows(t)] \nonumber \\ &= S(f) * F_s \sum_{m=-\infty}^\infty \delta(f - m F_s) \nonumber \\ &= F_s \sum_{m = -\infty}^{\infty} S(f - m F_s) \nonumber \end{align} \end{split}\]

Para entender la última igualdad revisemos la siguiente sección

¿Qué significa convolucionar con un tren de impulsos?¶

Matemáticamente la convolución entre dos señales discretas (de una dimensión) se define como

\[ (f * g)[n] = \sum_{m=-\infty}^\infty f[m] g[n-m] = \sum_{m=-\infty}^\infty f[n-m] g[m] \]

Podemos imaginar que $g$ se desplaza sobre $f$ (o viceverza)
En cada paso el $g$ desplazado se multiplica punto a punto con $f$ y luego se suman
El resultado de convolucionar $f$ con $g$ es una nueva función

Si $f$ es un tren de impulsos ocurrirá una «repetición» de $g$ (o viceverza)

La siguiente animación muestra en la imagen superior un pulso cuadrado que se desplaza sobre sren de impulsos en la figura superior. La figura inferior muestra el resultado de la convolución

%%capture
fig, ax = plt.subplots(2, figsize=(7, 4), sharex=True, tight_layout=True)
t = np.arange(-4, 4, step=1e-4)

def my_signal(t, a=0, T=0.5):
    s = np.zeros_like(t)
    s[np.absolute(t-a)<T] = 1
    return s

# 2 segundos de espacio entre los dientes de la peinte
def peineta(t):
    s = np.zeros_like(t)
    s[::20000] = 1 
    return s
    
conv_s = np.convolve(my_signal(t), peineta(t), mode='same')

def update(a = 0): 
    ax[0].cla(); ax[1].cla()
    p1, p2 = my_signal(t, 0.1*a - 4), peineta(t)
    ax[0].plot(t, p1, label='señal'); 
    ax[0].plot(t, p2, label='peineta'); 
    ax[0].legend()
    ax[1].plot(t, conv_s); 
    ax[1].set_xlabel('Tiempo [s]'); 
    ax[1].scatter(0.1*a -4, np.sum(p1*p2), s=100, c='k')
    return ()
    
anim = animation.FuncAnimation(fig, update, frames=80, interval=100, blit=True)

HTML(anim.to_html5_video())

Por lo tanto: Muestrear una señal en el tiempo hace que su espectro $S(f)$, sea cual sea, se repita infinitamente. Además el espectro se repite cada $F_s$ [Hz]

¿Cúal es el espectro de la señal Gaussiana?¶

La transformada de Fourier de

\[ s(t) = \exp \left ( -\frac{t^2}{2\sigma^2} \right) \]

es

\[ S(f) = \mathbb{FT}[s(t)] = \sqrt{2\pi}\sigma \exp \left ( -2\sigma^2 \pi^2 f^2 \right) \]

La función gaussiana es muy especial: El espectro de una gaussiana es otra gaussiana

Puedes ver la demostración de la transformada anterior aquí

En base a la gráfica siguiente notemos que

Mientras más «ancha» sea la gaussiana en el tiempo ($\sigma$ pequeño) más «angosto» será su espectro en frecuencia
Las gaussianas son del mismo ancho cuando $\sigma = \frac{1}{\sqrt{2\pi}}$

x = np.arange(-5, 5, step=0.001)
sigma = Slider(start=0.1, end=2, value=1/np.sqrt(2*np.pi), step=.01, title="sigma")

source = ColumnDataSource(data=dict(x=x, 
                                    gt=np.exp(-0.5*x**2/sigma.value**2),
                                    gf=np.exp(-2*(np.pi*x*sigma.value)**2)*np.sqrt(2*np.pi)*sigma.value
                                   ))

create_figure = lambda title : Figure(plot_width=350, plot_height=280, title=title, 
                                      toolbar_location="below", x_range=(-5, 5))
p1 = create_figure('Dominio de tiempo', )
p2 = create_figure('Dominio de frecuencia')
p1.line('x', 'gt', source=source, line_width=3)
p2.line('x', 'gf', source=source, line_width=3)

callback = CustomJS(args=dict(source=source, sigma=sigma), code="""
    var x = source.data['x'];
    var gt = source.data['gt'];
    var gf = source.data['gf']
    for (var i = 0; i < x.length; i++) {
        gt[i] = Math.exp(-0.5*Math.pow(x[i]/sigma.value, 2));
        gf[i] = Math.sqrt(Math.PI*2)*Math.exp(-2*Math.pow(Math.PI*x[i]*sigma.value, 2))*sigma.value;
        
    }
    source.change.emit();
""")


sigma.js_on_change('value', callback)
show(column(sigma, row(p1, p2)))

¿Cuál es el espectro de la gaussiana «discreta»?¶

Como ya vimos el espectro de una señal discreta es idéntico al de la señal continua pero «repetido» según el valor de la frecuencia de muestreo

Por ejemplo si $\sigma=1$ y $Fs = 2$ [Hz] tendríamos lo siguiente

def gaussiana(f, sigma):
    return np.sqrt(2*np.pi*sigma**2)*np.exp(-2*(np.pi*f*sigma)**2)

def espectro_discreto(f, sigma):
    S = np.zeros_like(f)
    for m in range(-20, 20):
        S += gaussiana(f - Fs*m, sigma)
    return S

def ventana(f, Fs):
    SW = np.zeros_like(f)
    SW[np.absolute(f) < Fs/2] = 1
    return SW

Fs = 2.
sigma = 1.
f = np.arange(-3*Fs, 3*Fs, step=1e-3)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma),  color='green', alpha=0.75,
        line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
        legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Ventana cuadrada')

p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

Por lo tanto: Si multiplicamos el espectro discreto (azul) por una ventana cuadrada de ancho $F_s$ (negro punteado) podemos recuperar el espectro original (verde) sin pérdidas

Esta es la base del siguiente teorema fundamental

Teorema del muestreo¶

Sea una señal continua $s(t)$ muestreada a $F_s$ [Hz] produciendo una señal digital $s[n] = s(t = n/F_s)$

Si

\[ f_{max} < \frac{F_s}{2}, \]

donde

$f_{max}$ es la componente de frecuencia más alta de la señal
$\frac{F_s}{2}$ es la frecuencia de Nyquist

entonces la señal analógica $s(t)$ puede ser recuperada perfectamente a partir de sus muestras discretas $s[n]$

Además el valor de la señal continua reconstruida es:

\[ s(t) = \sum_{n=-\infty}^{\infty} s[n] \text{sinc}(\pi F_s (t - n /F_s) ) \]

¿Qué pasa con el espectro «discreto» si no se cumple la condición anterior?¶

Asumamos que la frecuencia de muestre se mantiene $F_s=2$ [Hz] y que $\sigma$ disminuye

La frecuencia máxima de la gaussiana es la «última» frecuencia donde el espectro es distinto de cero

Si $\sigma$ disminuye la $f_{max}$ aumenta

Fs = 2.
sigma = 0.4
f = np.arange(-3*Fs, 3*Fs, step=1e-3)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma),  color='green', alpha=0.75,
        line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
        legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Ventana cuadrada')

p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

Si asumimos que $\sigma =1 $ se mantiene y que la frecuencia de muestreo disminuye se obtiene el mismo efecto

Fs = 0.75
sigma = 1.
f = np.arange(-3*Fs, 3*Fs, step=1e-3)

p1 = Figure(plot_width=600, plot_height=280, toolbar_location="below")
p1.line(f, gaussiana(f, sigma),  color='green', alpha=0.75,
        line_width=3, legend_label='Espectro original')
p1.line(f, espectro_discreto(f, sigma), line_width=3, alpha=0.75,
        legend_label='Espectro discreto')
p1.line(f, ventana(f, Fs), line_width=3, alpha=0.75, color='black',
        line_dash='dashed', legend_label='Ventana cuadrada')

p1.xaxis[0].axis_label = 'Frecuencia [Hz]'
show(p1)

Debido al solapamiento en el espectro discreto (azul) se vuelve imposible recuperar el espectro original (verde) sin alteraciones. Podríamos recuperar su parte central, pero eso significa perder información de alta frecuencia

El solapamiento espectral se llama aliasing

Aliasing¶

El espectro de una señal muestreada es periódico en $F_s$
Si originalmente la señal tenía componentes con frecuencias mayores a $\frac{F_s}{2}$ se produce un «traslape» o «solapamiento» espectral
Este fenomeno se llama aliasing y los componentes traslapados se denominan aliases

Cuando existe aliasing no se puede reconstruir la señal original

¿Cómo eliminamos el aliasing?¶

Necesitamos que todas las frecuencias del espectro sean menores que $\frac{Fs}{2}$

Podemos lograr esto mediante

Filtrado: Eliminar las frecuencias mayores a $\frac{F_s}{2}$ (Próxima unidad)
Aumentar $F_s$ tal que sea dos veces mayor que la frecuencia máxima de interés, no siempre es fácil saber cuál será su valor a priori

Ejemplo¶

Se tiene una señal con tres componentes a 1Hz, 12Hz y 23Hz, respectivamente
Modifica la frecuencia de muestreo y detecta los componentes erroneos o aliases en el espectro
Verifica como afectan a la reconstrucción de la señal

import scipy.fft as fftpack

fig, ax = plt.subplots(4, figsize=(6, 6));
t = np.arange(-3, 3, step=0.01); N = len(t)
s =  np.cos(2*np.pi*t) - 0.5*np.sin(2.0*np.pi*12*t + np.pi/3) + 0.25*np.sin(2.0*np.pi*23*t) 
ax[0].plot(t, s, '-'); ax[0].set_title('Original')

def update(Fs):    
    t_short = t[::int(100/Fs)]; s_short = s[::int(100/Fs)]
    ax[1].cla(); ax[1].plot(t_short, s_short, '.'); ax[1].set_title('Sampled at %0.2f' %(Fs))
    S = fftpack.fft(s_short); f = fftpack.fftshift(fftpack.fftfreq(n=len(S), d=1.0/Fs)); 
    ax[2].cla(); ax[2].plot(f, fftpack.fftshift((np.absolute(S))), linewidth=2);
    S_pad = np.concatenate((S[:int(len(S)/2)], np.zeros(shape=(N-len(S))), S[int(len(S)/2):]))
    ax[2].set_title('Spectrum '); s_recon = np.real(fftpack.ifft(S_pad))
    ax[3].cla(); ax[3].plot(t, s_recon*len(t)/len(t_short), '-', linewidth=2); 
    ax[3].set_title("Reconstruction from padded spectrum");
interact(update, Fs=SelectionSlider(options=[5, 10, 25, 33, 50, 100], value=100));

cosas viejas

Principio o Teorema de incertidumbre¶

El principio de incertidumbre de Heisenberg nos dice que la precisión (certeza) con que medimos la posición de una particula es inversamente proporcional a la precisión con que medimos su momentum lineal:

\[ \Delta x \Delta p \geq \frac{h}{4\pi}, \]

donde $h$ es la constante de Planck.

En señales existe un principio análogo: No podemos especificar con infinita precisión la localización temporal y frecuencial de una señal al mismo tiempo.

Ejemplo: Pulso cuadrado¶

Teníamos

\[\begin{split} s(t) = \begin{cases} 1, & |t| < T \\ 0, & |t| > T\end{cases} \end{split}\]

con

\[ S(\omega) = \mathbb{FT}[s(t)] = 2T \text{sinc}(\omega T) \]

y notamos que a mayor ancho en tiempo menor ancho en frecuencia y viceversa.

Definiendo la resolución o ancho de la señal como la distancia entre sus cruces por cero, es fácil verificar que:

$\Delta t = 2T$
$\Delta \omega = \omega|_{\text{sinc}(- \omega T) = 0}^{\text{sinc}(\omega T) = 0} = \frac{2\pi}{T}$
$\Delta t \Delta \omega = 4 \pi$

Es decir a menor ancho en el tiempo, mayor ancho en frecuencia y viceverza

plt.close('all'); fig, ax = plt.subplots(1, 2, figsize=(6, 3))
t = np.linspace(-5, 5, num=500); f = np.linspace(-5, 5, num=500)
ax[0].set_xlabel('Time [s]'); ax[1].set_xlabel('Frequency [Hz]');
line_square = ax[0].plot(t, np.zeros_like(t), linewidth=4); ax[0].set_ylim([-.1, 1.1])
line_gabor1 = ax[0].plot([-1, 1], [0, 0], 'r--')
line_sinc = ax[1].plot(f, np.zeros_like(f), linewidth=4)
line_gabor2 = ax[1].plot([-0.5, 0.5], [0, 0], 'r--')

def update(T):    
    s = np.zeros_like(t); s[(t> - T) & (t<  T)] = 1    
    line_square[0].set_ydata(s);  line_gabor1[0].set_xdata([-T, T])
    S = 2*T*np.sinc(2*f*T); line_sinc[0].set_ydata(S); 
    line_gabor2[0].set_xdata([-0.5/T, 0.5/T])
    ax[1].set_ylim([-2*T/np.pi, 2.2*T])
    
interact(update, T=SelectionSlider(options=[1/8, 1/4, 1/2, 1, 2, 4], value=1));

Intuición: Resolución frecuencial de una sinusoide¶

from matplotlib import patches

plt.close('all'); fig, ax = plt.subplots(figsize=(7, 3))
t = np.arange(-4, 4, step=1e-2)
def update(dt):
    ax.cla(); ax.plot(t, np.cos(2.0*np.pi*t))
    ax.plot(t, np.cos(2.0*np.pi*(1.0+0.125/dt)*t))
    rect = patches.Rectangle((-dt, -1), 2*dt, 2, fill=False, lw=2)
    ax.add_patch(rect)
interact(update, dt=SelectionSlider(options=[0.5, 1., 2.]));

Denis Gabor (1946) fue el primero en darse cuenta de que el principio de incertidumbre aplica para señales. Formalmente su teorema:

Para una señal con energía finita $$ E = \int |s(t)|^2 dt $$ con valor medio $$ \langle t \rangle = \frac{1}{E} \int t |s(t)|^2 dt, $$ y varianza temporal $$ (\Delta t)^2 = \frac{1}{E} \int (t - \langle t \rangle)^2 |s(t)|^2 dt, $$ cuya transformada de Fourier $\mathbb{FT}[s(t)] = S(\omega)$ tiene un valor medio en frecuencia $$ \langle \omega \rangle = \frac{1}{E} \int (\omega - \langle \omega \rangle) |S (\omega)|^2 d \omega $$ y varianza frecuencial $$ (\Delta \omega)^2 = \frac{1}{E} \int (\omega - \langle \omega \rangle)^2 |S(\omega)|^2 d\omega $$

Entonces se cumple que

\[ \Delta t \Delta \omega \geq \frac{1}{2}, \]

es decir $\Delta t$ y $\Delta \omega$ no pueden ser arbitrariamente pequeños. Esto se conoce también como límite de Gabor

Ejemplo: Transformada de Fourier de una Gaussiana¶

La función Gaussiana se define como

\[ s(t) = \alpha e^{- \beta t^2} \]

y su transformada de Fourier es

\[\begin{split} \begin{align} S(\omega) &= \alpha \int e^{- \beta t^2} e^{-j\omega t} dt \nonumber \\ &= \alpha e^{- \frac{\omega^2}{4\beta}}\int e^{-\beta (t + \frac{j\omega}{2\beta})^2 } dt \nonumber \\ &= \alpha e^{- \frac{\omega^2}{4\beta}} \sqrt{\frac{\pi}{\beta}} = \hat \alpha e^{- \hat \beta \omega^2}, \nonumber \end{align} \end{split}\]

es decir otra gaussiana.

La gaussiana es la unica función que cumple $\Delta t \Delta \omega = \frac{1}{2}$

Definiciones¶

Función limitada en el tiempo $$ s(t) = 0 \quad \forall |t| > T, $$ para alguna constante $T$
Función limitada en ancho de banda $$ S(\omega) = 0 \quad \forall |\omega| > \Omega, $$ para alguna constante $\Omega$

Enventanado¶

En la práctica no trabajamos con señales de duración infinita
Una señal de duración infinita puede hacerse finita multiplicando por una ventana finita
Por ejemplo $$ s_T(t) = \cos(\omega_0 t) \text{rect}(t/T), $$ donde $$ \text{rect}(x) = \begin{cases} 1 & |x| \leq 1 \\ 0 & |x| > 0 \end{cases} $$ es una ventana rectangular, i.e. pulso cuadrado

¿Cúal es la transformada de Fourier de $s_T(t)$?

\[\begin{split} \begin{align} S_T(\omega) &= \int s_T(t) e^{-j\omega t} dt \nonumber \\ &= \int_{-T}^T \cos(\omega_0 t) \cos(\omega t) dt \nonumber \\ &= T \text{sinc}((\omega - \omega_0)T) + T \text{sinc}((\omega + \omega_0)T) \nonumber \end{align} \end{split}\]

Donde usamos que $$ 2 \cos(\alpha)\cos(\beta) = \cos(\alpha - \beta) + \cos(\alpha + \beta) $$

Que también se resuelve usando la propiedad de modulación $$

()¶\[\begin{align} S_T(\omega) &= T \text{sinc}((\omega - \omega_0)T) + T \text{sinc}((\omega + \omega_0)T) \nonumber \\ &= T \text{sinc}(\omega T) * \left[ \delta(\omega - \omega_0) + \delta(\omega + \omega_0) \right] \nonumber \\ &= \frac{1}{2\pi} \mathbb{FT}[\text{rect}(t/T)] * \mathbb{FT}[\cos(\omega_0 t)] \nonumber \end{align}\]

$$

Es decir que multiplicar por una ventana rectangular en el tiempo es equivalente a convolucionar con un sinc en frecuencia
¿Qué repercusión tiene esto en el espectro?

plt.close('all'); fig, ax = plt.subplots(2, figsize=(6, 4), tight_layout=True)
f = np.linspace(-3, 3, num=500)
def update(T, window):
    t = np.arange(-T, T, step=1e-2); 
    tf = np.cos(2.0*np.pi*t[:, np.newaxis]*f[:, np.newaxis].T);
    if window == "rect":
        w = np.ones_like(t)
    elif window == "gaussian":
        w = np.exp(-0.5*t**2/(0.333*T)**2)
    s = w*np.cos(2.0*np.pi*t); 
    S = np.average(s[:, np.newaxis]*tf, axis=0)
    ax[0].cla(); ax[1].cla(); ax[0].set_title(r'$s_T(t) = rect(t/T)cos(2\pi t)$')
    ax[0].plot(t, s); ax[1].plot(f, S); ax[1].set_title(r'$\Re[S_T(f)]$')
interact(update, T=SelectionSlider(options=[0.5, 1, 2, 4, 8, 16], value=8),
         window=SelectionSlider(options=["rect", "gaussian"]));

Fuga espectral¶

La fuga espectral (spectral leak):

En una distribución de la energía de un cierto componente espectral hacia frecuencias vecinas
Aparece como lobulos laterales en torno a los componentes espectrales
Sólo aparece si la señal es no periódica dentro del rango de observación
Está asociado a las discontinuidades provocadas por los bordes de la señal
- Una discontinuidad fuerte puede ocupar mucho espacio en frecuencia ¿Por qué?

Multiplicando la señal por una ventana apropiada podemos reducir este efecto a un rango más acotado de frecuencia

Ventana de Hann¶

Se define como $$ w[n] = 0.5 - 0.5 \cos \left( \frac{2\pi n}{N-1} \right), \quad n=0, 1, \ldots, N-1, $$ y tiene bordes que tienden a cero

plt.close('all'); fig, ax = plt.subplots(2, figsize=(6, 4))
def update(f0, window, log_scale=False):
    n = np.arange(0, 100);
    if window == "rect":
        w = np.ones_like(n)    
    elif window == "hamming":
        w = 0.5 - 0.5*np.cos(2*np.pi*n/(len(n)))
    s = w*np.cos(2.0*np.pi*f0*n); f = fftpack.fftshift(fftpack.fftfreq(n=len(s), d=1))
    S = fftpack.fftshift(np.absolute(fftpack.fft(s)))
    ax[0].cla(); ax[1].cla(); ax[0].set_title(r'$w[n]cos(2\pi f_0 n)$')
    ax[0].plot(n, s); ax[1].plot(f, S); ax[1].set_title(r'$|S_T(f)|$')
    if log_scale:
        ax[1].set_yscale('log')
interact(update, f0=SelectionSlider(options=[0.05, 0.05214178, 0.055414684, 0.2, 0.241245], value=0.05),
         window=SelectionSlider(options=["rect", "hamming"]));

Toma 2: Efectos de la frecuencia de muestreo y de la ventana de observación en el espectro¶

Fila superior, linea azul: señal con frecuencia de muestreo 50 Hz y duración 10 segundos
Superior izquierda (cruces rojas): señal con un décimo de la frecuencia de muestreo pero igual largo temporal
Superior derecha (cruces rojas): señal con igual frecuencia de muestreo pero con un décimo del largo temporal

La frecuencia de muestreo influye en la frecuencia máxima que podemos estudiar (frecuencia de Nyquist)

El largo temporal o ventana de observación influye en la resolución frecuencial del espectro

f0=1.234
t_plot = np.linspace(0, 10, 1000)

x_plot = np.cos(2.0*np.pi*f0*t_plot) + 0.4*np.cos(2.0*np.pi*2*f0*t_plot +np.pi/4) + 0.1*np.random.randn(1000)

t, x = t_plot[::10], x_plot[::10]
X = fftpack.fft(x)[:len(x)//2]
freq = fftpack.fftfreq(n=len(x), d=t[1]-t[0])[:len(x)//2]

fig, ax = plt.subplots(2, 2, figsize=(7, 4), tight_layout=True)
ax[0, 0].plot(t_plot,x_plot)
ax[0, 0].scatter(t, x, marker='x', c='r')
ax[1, 0].plot(freq, np.absolute(X))

t, x = t_plot[:100], x_plot[:100]
X = fftpack.fft(x)[:len(x)//2]
freq = fftpack.fftfreq(n=len(x), d=t[1]-t[0])[:len(x)//2]
ax[0, 1].plot(t_plot,x_plot)
ax[0, 1].scatter(t, x, marker='x', c='r')
ax[1, 1].plot(freq, np.absolute(X));

Toma 2: Incertidumbre frecuencia-tiempo¶

Mientras más pequeño es la «separación en frecuencia» de dos señales, más tiempo debemos observarlas para poder diferenciarlas

t, dt = np.linspace(0.0, 5.0, num=1000, retstep=True)

fig, ax = plt.subplots(4, figsize=(8, 6), tight_layout=True)

for k in range(4):
    ax[3-k].plot(t, np.cos(2.0*np.pi*2.1*t))
    ax[3-k].plot(t, np.cos(2.0*np.pi*2.1*(1+0.01+0.01*k)*t))
    ax[3-k].set_title('Error relativo: %0.2f' %(0.01+0.01*k))

INFO183 Análisis de sistemas lineales

Efectos del muestreo y aliasing¶